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Abstract. 

We present a numerical method, based on a FEM simulation, for the determination 
of the gravitational field generated by massive objects, whatever geometry and space 
mass density they have. The method was applied for the determination of the self 
gravity effect of an absolute cold atom gravimeter which aims at a relative uncertainty 
of 10~^. The deduced bias, calculated with a perturbative treatment, is finally 
presented. The perturbation reaches (1.3 ±0.1) x 10~^ of the Earth's gravitational 
field. 
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1. Introduction 

Gravitational forces from surrounding masses are usually negligible compared to electro- 
magnetic forces and thus not considered as relevant in laboratory based experiments. 
Nevertheless, they can significantly impact on accurate forces measurement [Ij or even 
constitute the dominant effect [2, 3j. As an example, an aluminum disk with 50 cm 
diameter and 4 cm thickness generates a gravitational field ranging from 4 to 1 parts in 
10^ of the Earth's gravity field from its surface to 30 cm along its axis. Components 
with similar mass density and dimensions are commonly used in precise instruments 
developed for measuring the gravity field with uncertainties of few parts in 10^. The 
measurement principle of these instruments is based on the determination of the ballistic 
trajectory followed by a mass, i.e. corner cube retro-refiectors and recently atoms, during 
their vertical free-fall under the infiuence of gravity [H El Ej . Due to constrains in their 
design, their mass distribution is in general not symmetric with respect to the trajectory 
of the test mass. The effect of parasitic attractions which is called self gravity effect 
when restricted to the infiuence of the device itself, has to be carefully estimated and if 
necessary corrected for in order to guarantee the accuracy of the measurement. 

In the case of the absolute corner cube gravimeter FG5, authors in [4J indicate 
that the g measurement is corrected from the attraction of the apparatus, with an 
uncertainty of 10~^ m.s~^ without giving details on its calculation. A similar level of 
uncertainty is also necessary in the case of the cold atom gravimeter (GAG) described 
in [7j , currently operating and improved within the framework of the LNE watt balance 
experiment [HJ [9] . 

In this paper we describe a numerical method for computing the gravitational field 
generated by extended and continuous massive objects, whatever geometry and space 
mass density they have. Application of the method for the quantitative estimate of the 
self gravity error on the GAG is also given. 

2. Analogy between gravitational and electrical interactions 

Mass and electric charge are respectively the sources of the gravitation and 
electromagnetic interaction between bodies jlO]. In particular the Newton's and 
Goulomb's laws which describe respectively the gravitational and electric force occurring 
between two masses and two charges, have the same behavior. This analogy can be 
exploited to compute gravitational forces using methods originally developed for electro- 
magnetic forces, by replacing individual charges by individual masses, space charge 
density by space mass density and I/Attcq by G where eo is the vacuum permittivity and 
G is the gravitational constant. 

In general, the gravitational field can be computed by solving cumbersome 
integrals of three-dimensional vector functions. Analytical solutions exist only for 
axial-symmetric and homogeneous bodies, e.g. a spherical shell, a solid sphere, a right 
rectangular prism, a right polygonal prism and a polyhedron [11]. For more complex 
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systems, which can not be easily broken down in simphfied parts, the Finite Element 
Method (FEM) helps in finding an approximate numerical solution of the integrals. The 
essential characteristic of this method is the mesh discretization of a continuous domain 
in a set of discrete sub-domains. For this application we used the FEM software packages 
developed by COMSOL Multiphysics [12], in particular the electro-static module. 

3. Finite element method and expected accuracy 

In our application, the three-dimensional model of the object in a Oxyz frame of 
reference defines a body domain BD fiUed-in a given mass density, whereas the 
surrounding domain SD shapes a vacuum space where we are interested in calculating 
the gravitational field. The interface surface Sj between BD and SD corresponds to 
the body geometry while the external surface Se limiting SD depends on the boundary 
conditions. 

Assigning the equipotential condition to Se simplifies the problem. Indeed, 
equipotential surfaces are strongly dependent on the body geometry but at infinite 
distance they become spheres centered on the center of mass CM (see Figure [l] 
for an example concerning an ellipsoidal body). Knowledge about the geometry of 
equipotential surfaces near complex objects is usually poor and the extension of SD 
to infinity is not numerically possible. Nevertheless, defining Se by a sphere centered 
on CM is the best choice, independently on the object. The systematic error arising 
from the possible field distortion can be minimized by increasing the sphere radius. The 
upper limit in stepping up the radius of Se is set by the finite number of mesh elements 
allowed by the memory of the solver. 

In fact, increasing the volume of SD results in a higher scattering of the solution due 
to a lower discretization of the space where the gravity field is computed. In conclusion 
the final uncertainty depends on a compromise between the bias due to the boundary 
condition and the scattering due to the mesh size. 

COMSOL Multiphysics offers several ways to set the mesh size. In case of objects 
having hollow or highly sharp geometries the algorithm must resolve the gravitational 
filed in great detail only on some portions of the domain. In this situation best results 
can be achieved using the adaptive mesh generation. After a first solution obtained 
with the mesh size selected by the user, the algorithm provides a second solution after 
optimizing the mesh size in those regions requiring an higher resolution. Moreover, the 
ultimate accuracy compatible with the memory of the solver is achieved by constraining 
an extra-fine mesh in those regions where the solution is required. 

To evaluate the expected uncertainty of the FEM method, a comparison was carried 
out between the analytical solution and the numerical results of a simulation of the same 
problem. Lets consider a right rectangular prism with uniform density and centered at 
the origin of a Oxyz frame of reference, with the z axis oriented upwards. Sides of the 
prism are parallel to x, z axes with dimensions equal to 2L, 2L and L, respectively. 
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Figure 1. Equipotential surfaces generated by an ellipsoidal body in its surrounding. 
At sufficiently large distances from the object, these equipotential are spherical. 

The z component of the gravitational field generated at a point (x, z) is pTl [T3] : 

2 2 2 / X 

r = -Gpm XI 5Z 5Z 1^'^^ ( ^^^y^ + ^'^^"^ + ^^^^^ + ^'^^"^ ~ ^^^^^^ ) (i 



where Xi x - L, yi y - L, Zi z - L/2, rijk ^x'j + y'j + zl and 

l^ijk = (— 1)^. Figure [2] represents the expected gravitational field along 
two lines parallel to the z axis, within a range equal to ±20L. The black and gray 
curves correspond to two vertical lines passing through points (0,0,0) and (2L,0,0), 
respectively. Data are normalized with the maximum vertical field T^ax^ occurring at 
point (0,0, —L/2). The FEM simulation was performed by limiting SD with a sphere 
having a radius equal to twenty times L. Both the domains were meshed with an 
adaptive mesh generation after uniform constrain of 15 mesh elements every L on the two 
lines where the solution is needed. The relative difference between the simulation results 
and the analytical solution are shown in Figure [3| As in the previous figure, the black 
curve concerns the line passing through (0,0,0) and the gray curve concerns the line 
passing through (2L,0,0). The systematic error occurring near Se and the scattering 
near the body geometry are kept below 0.1 % of the generated field. Significantly 
better results can be achieved by running the FEM algorithm on a platform with better 
memory power, which allows to increase the radius of Se and the number of mesh 
elements constrained along the lines. 



Perturbations of the local gravity field due to mass distribution on precise measuring instruments: 




Figure 2. Gravitational field of a prism of dimensions 2L, 2L, L, calculated with the 
FEM simulation, along two lines parallel to the z axis passing at the center of the 
prism (black line) and at the distance L of the prism (gray line). 
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Figure 3. Relative difference results between the FEM simulation of a prism 
(2L, 2L, L) represented in Figure [2j and the analytical solution according to equation [l] 
for the two vertical lines passing at the center of the prism (in black) and L of the 
prism (in gray). 



4. The cold atomic gravimeter 

A scheme of the C AG setup is presented in [Gj [7] . It performs a cyclic measurement of 
the gravity acceleration g with a cloud of ^''Rb cold atoms used as a test mass [M]. The 
gravimeter core is shown in Figure [4} 

A detailed description of the principle of the gravimeter can be found in [14j . Briefly, 
an atomic cloud is loaded in a 3D-M0T (Magneto Optical Trap) and is further cooled 
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down to 2 /xK before being released. While the atoms fall down, three Raman pulses 
separated by T (7r/2 — tt — 7r/2) split, redirect and recombine the atomic wave packets. 
They are induced by two vertical counterpropagating laser beams of wave-vectors ki 
and k2 which couple the hyperfine levels |F = 1) and |F = 2) of the ^5^1/2 ground state 
via a two-photon transition [15j. They are delivered to the atoms through a single 
collimator [14j and retro-reflected with a mirror placed inside the vacuum chamber. 
Due to conservation of angular momentum and to the Doppler shift induced by the 
free fall of the atoms, only two counter-propagating of the four beams will drive the 
Raman transitions. This feature allows to perform interferometers using effective wave- 
vector fceff = fci — ^2 pointing upwards or downwards. Finally, thanks to the state 
labelling method [16j, the interferometer phase shift A$ which is the diflFerence of the 
atomic phases accumulated along the two paths / and // (Figure [i]), is deduced from 
a fluorescence measurement of the populations of each of the two states. It is given 
by [E!: 

A$ = ±|4ffll^|r2 (2) 

where |/Ceff| = + |/c2| for counter-propagating beams. Performing the interferometer 
with initial atoms in state |F = 1) leads to different paths / and // using fceff^ or fceff^ 
(in black and gray on Figure [4]). The interferometer takes place in between the center 
of the MOT chamber (z = m) and the detection setup area {z = —0.16 m), through 
the free-fall vacuum chamber (FF). The atoms travel in an empty vertical cylinder of 
40 mm diameter. The actual trajectory depends on the delay ti of the flrst pulse and 
on the time T separating the Raman pulses. In our case, ti = 16 ms and T = 70 ms, so 
that the interferometer occurs in the region [-1.3 mm; -120.2 mm]. 



5. Mass attraction along the trajectory of the gravimeter 

The GAG apparatus consists of four main parts, (i) the gravimeter core presented 
in Figure [ij (n) an isolation platform used to fllter the vibrations due to background 
noise [18j, (Hi) a thermal and acoustic insulation box and (iv) an optical bench [20j . 

The level of details required in the modeling of each of the four parts depends on 
their mass and location with respect to the atomic trajectories. Most of the attention 
was focused on the gravimeter core, made of several sub-components located close to 
the free-falling atoms. The isolation platform was modeled with an homogeneous box 
approximating the inner mass distribution. Due to the distance to the atoms, details 
of the elements inside the platform would not change the results. The insulating box 
was easily modeled with four lateral walls and a roof. Although the mass of the optical 
bench is signiflcant, it is far away and the direction of its gravitational attraction is 
nearly perpendicular to the free falling motion. The numerical simulation of its effect 
was not necessary (maximum effect < lO""*^^ m.s~^). 

Taking advantage of the linear additivity of the gravitational fleld, the perturbation 
was evaluated by splitting the whole setup on its sub-components, computing the 
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Figure 4. Scheme of the interferometer. Left, vacuum chamber partly cut together 
with some parts of the apparatus, with the magnetic shields partially removed. Center, 
enlarged view of the free fall section. Right, atomic trajectories with Raman wave- 
vector pointing either downwards (black line) or upwards (gray lines). The Raman 
pulses are depicted with gray filled areas. For the sake of clarity the value of the 
wave- vector /ceff has been multiplied by 15 to distinguish the separation between path 
/ (OABD) (dash lines) and path // (OACD) (continuous lines). 

vertical components of the field intensities along the atoms trajectory and adding the 
contributions. 

All the objects and the results were modeled in a Oxyz reference frame, having 
the origin located at the start of the atoms drop, with the z axis vertical and upwards 
oriented. 

We choose to use a Se radius 100 times times larger than the maximum distance 
between the object surface and its center of mass, significantly larger than the factor 20 
used in section |3] for which the maximum error was as low as 0.1%. On the other hand, 
the adaptive mesh and the elements constrain along the trajectory were used depending 
on the magnitude of the effect. In most cases, user selection of the mesh size without 
element constraint were enough, from numerous tests carried out for the evaluation of 
the modeling uncertainty we estimate the net error well below 1 % of the generated 
gravity field. 

Figure [5] shows the average effect for each components, from the biggest to the 
lowest one. The larger attraction is due to the free fall chamber which increases the 
g value by more than 3 x 10~^ m.s~^. The magneto optical trap (MOT) vacuum 
chamber compensates for most of this effect as anticipated during the design. The 
third component modifying g with an impact larger than 10~^ m.s~^ is the detection 
setup, other elements having a smaller impact. As an example, the insulating box acts 
upwards by about 0.3 x 10~^ m.s~^. Figure [6] displays the perturbation due to the mass 
attraction along the trajectory of the atoms. The parasitic acceleration is downwards 
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Figure 5. Average gravity effect for sub- components of the CAG. 




oriented and varies from 1.3 x 10~^ m.s~^ to 2.9 x 10~^ m.s~^ along the atoms trajectory 
with a minimum of 0.5 x 10~^ m.s~^ between the second and the third Raman pulse. 
The effect is stronger after the last pulse, in the detection area (< — 12 mm), and close 
to 3 X 10"^ m.s"^ 

6. Effect on the gravity measurement 

When small enough, parasitic phase shifts of the interferometer can be accurately 
calculated using a perturbative path integral treatment [21j. This approach already 
used in [22l [23] to treat the effect of a linear vertical gravity gradient, is used here to 
calculate the perturbation of the interferometer phase due to the mass attraction effect 
of the device itself (A$r)- It consists on integrating the perturbed Lagrangian Cpert 
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along the unperturbed paths / and //: 



[z{t),z{t)]dt 



(3) 



where h is the reduced Planck constant. In our case, the perturbed Lagrangian is hnked 
to the perturbed potential energy which is obtained by integrating the attraction force: 



>Cpert = -^Ep = ^ J r(z)(iz mf{z) 



(4) 



where m is the mass of the falling atoms. Figure [7] shows the function f{z) obtained 
integrating the mass attraction T{z) plotted in Figure [6} 

The position of the falling atoms z{t) is calculated for each path / and // taking 
into account the Raman pulses changing the velocity in B for path / and in A and C 
for path //, leading to the functions represented on Figure |4| 



ZAB{t), ZBD{t), ZAc{t) and ZcD{t) 
Combining equations |4] and [5] in equation [3| A$r can be expressed as 



(5) 



m 



/ f{zAB{t))dt- / f{zBD{t))dt 

\Jti Jti+T 
ti+T pti+2T 

f{zAc{t))dt+ / f{zcD{t))dt 



with ti the time of the first 
bias obtained with equation 



Raman pulse. With our experimental parameters, the 
2] is 1.27 X 10~^ m.s~^ whatever the direction of fcefr- The 
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difference between the two directions is negligible (10~^^ m.s~^). Moreover, due to 
the finite temperature, the atoms also move radially while falling. At 2 /xK, a point- 
like atomic cloud reaches a diameter of 4 mm in the detection zone. We thus 
calculate the mass attraction F for trajectories off-centered by 5 mm in the two horizontal 
directions x and ^, to calculate the corresponding shifts on the g measurement. We found 
differences with respect to the centered trajectory, as small as 4 x 10~^^ m.s~^. 

The global uncertainty in the calculation can be estimated by summing 
quadratically the uncertainties corresponding of the individual pieces. We find 5 x 
2Q-10 g-2 addition to the uncertainty in the calculation, which can be further 
reduced improving the mesh, we have also to account for additional uncertainties due 
to modeling in the mass distribution. They arise from simplifications in the geometrical 
description of the considered sub-components, from the infiuence of neglected sub- 
components and from imperfect knowledge of the densities. We assign a relative 
uncertainty of 1% as a conservative estimate for these contributions which is still 
negligible. The total uncertainty for the self gravity error is therefore increased to 
10~^ m.s~^ which is, up to now, much lower than the targeted accuracy of the CAG. 

Finally we consider the g data collected with the CAG to be corrected from 
(1.3 ± 0.1) X 10~^ m.s~^ due to the self gravity effect. This result agrees with the 
preliminary rough estimation of the effect of (0 ± 2) x 10~^ m.s~^ performed for last 
ICAG'09 [24]. 

7. Conclusion 

We have described a numerical method based on a FEM simulation for accurate 
determination of the gravitational attraction generated by any distribution of massive 
objects, applied to the self gravity effect on an atomic gravimeter. Using a perturbative 
treatment we find a relatively small systematic effect of (1.3 ±0.1) x 10~^ m.s~^ thanks 
to a symmetric design of the vacuum chamber. It shows this study was required in order 
to achieve the targeted accuracy. 

This method allows for an accurate determination of the attraction effect of any 
distribution of surrounding masses with an uncertainty arbitrary low providing densities 
and geometries are perfectly known. It can be applied to a large class of experiments 
such as watt balances [25l [26l [27j , torsion balances j2l EHJ [29l [30] or atom interferometers 
for G measurements [3] and experiments testing the equivalence principle [3l] , for which 
very accurate determination of the parasitic forces are mandatory in order to reach the 
ultimate level of performance. 
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